##### SLUMP AND EXPANSION PROBABILITY FROM AUERBACH & GORODNICHENKO (2013) ####

# for each country: 

# F(z_{t}) = exp(-\gamma*z_{t}) / (1 + exp(-\gamma*z_{t})), \gamma > 0,
# var(z_{t}) = 1, E(z_{t}) = 0


macro_identify <- macro_data

state_identifier <- function(df, var, countryname, thresh, tol){
  
  # select and filter dataset 
  
  df <- as.data.frame(df)
  
  df <- df %>% dplyr::filter(country==countryname) %>%
    dplyr::select(year,var)
  
  df <- df[complete.cases(df),]
  
  # make variable scale invariant with z-scores
  
  df$mean_var <- mean(df[,2])
  df$sd_var <- sd(df[,2])
  
  df$std_var <- (df[,var] - df$mean_var)/df$sd_var
  
  # set gamma probability space and run through dataset with logistic function
  # I generate 1000 random numbers between 0.1 and 5 to minimize gamma
  # to fulfill the condition set by the threshold
  
  gamma <- runif(1000, 0.1, 5) 
  
  hvec <- as.data.frame(sapply(1:1000, function(x)exp(-gamma[x]*df$std_var) / (1 + exp(-gamma[x]*df$std_var))))
  
  # trim observations at threshold level into binary indicators
  
  df[,6:1005] <- as.data.frame(sapply(1:1000, function(i)ifelse(hvec[,i] > thresh, 1, 0)))
  
  # clean dataset for merging
  
  df <- df %>% dplyr::select(year,starts_with("V"))
  
  df <- reshape2::melt(df, id="year")
  
  df <- df %>% group_by(variable) %>%
    dplyr::mutate(value_mean=mean(value)) %>%
    ungroup()
  
  # filter columns which fulfill threshold condition for state distinction
  
  thresh_inv = 1-thresh
  
  df <- df[df$value_mean>=(thresh_inv-tol) & df$value_mean<=(thresh_inv+tol),]
  
  df <- df %>% group_by(year) %>% dplyr::filter(value_mean==min(value_mean)) %>%
    dplyr::distinct(year,value_mean,.keep_all=T) %>%
    dplyr::select(year,variable,value) %>% rename(state=value)
  
  # prepare merging with weights on state
  
  hvec$year <- df$year
  
  hvec <- reshape2::melt(hvec, id="year")
  
  df <- left_join(df, hvec, by=c("year","variable"))
  
  df$country <- countryname
  
  df <- df %>% dplyr::select(!variable) %>%
    rename(weight=value) %>% relocate(country,.before=year)
  
  return(df)
}

## Following Auerbach and Gorodnichenko, I manually adjust the parameter gamma so every country
## is 20% of the time in a deficit episode; of course, the assumption of equal deficit episodes
## is not borne out by the facts; some countries ran deficits with a much higher frequency than
## others. However, this method still identifies the years where a deficit regime has been 
## *most likely*. 


## Real GDP Growth to identify slumps ##

au1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Australia", thresh=0.8, tol=0.01)
at1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Austria", thresh=0.8, tol=0.01)
be1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Belgium", thresh=0.8, tol=0.01)
ca1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Canada", thresh=0.8, tol=0.01)
dk1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Denmark", thresh=0.8, tol=0.01)
fi1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Finland", thresh=0.8, tol=0.01)
fr1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="France", thresh=0.8, tol=0.01)
de1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Germany", thresh=0.8, tol=0.01)
ie1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Ireland", thresh=0.8, tol=0.01)
it1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Italy", thresh=0.8, tol=0.01)
jp1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Japan", thresh=0.8, tol=0.01)
pt1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Portugal", thresh=0.8, tol=0.01)
es1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Spain", thresh=0.8, tol=0.01)
se1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Sweden", thresh=0.8, tol=0.01)
uk1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="United Kingdom", thresh=0.8, tol=0.01)
us1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="USA", thresh=0.8, tol=0.01)

recdata1 <- as.data.frame(rbind(au1,at1,be1,ca1,dk1,fi1,fr1,de1,ie1,it1,jp1,pt1,es1,se1,uk1,us1))

ggplot(data=recdata1, aes(x=year, y=weight)) + geom_line() + facet_wrap(facets="country")

recdata1 <- recdata1 %>% dplyr::rename(recession=state, recession_weight=weight)

macro_data <- left_join(macro_data, recdata1, by=c("country","year"))


## Deficits ##

au1 <- state_identifier(df=macro_identify, var="deficit",countryname="Australia", thresh=0.8, tol=0.01)
at1 <- state_identifier(df=macro_identify, var="deficit",countryname="Austria", thresh=0.8, tol=0.01)
be1 <- state_identifier(df=macro_identify, var="deficit",countryname="Belgium", thresh=0.8, tol=0.01)
ca1 <- state_identifier(df=macro_identify, var="deficit",countryname="Canada", thresh=0.8, tol=0.01)
dk1 <- state_identifier(df=macro_identify, var="deficit",countryname="Denmark", thresh=0.8, tol=0.01)
fi1 <- state_identifier(df=macro_identify, var="deficit",countryname="Finland", thresh=0.8, tol=0.01)
fr1 <- state_identifier(df=macro_identify, var="deficit",countryname="France", thresh=0.8, tol=0.01)
de1 <- state_identifier(df=macro_identify, var="deficit",countryname="Germany", thresh=0.8, tol=0.01)
ie1 <- state_identifier(df=macro_identify, var="deficit",countryname="Ireland", thresh=0.8, tol=0.01)
it1 <- state_identifier(df=macro_identify, var="deficit",countryname="Italy", thresh=0.8, tol=0.01)
jp1 <- state_identifier(df=macro_identify, var="deficit",countryname="Japan", thresh=0.8, tol=0.01)
pt1 <- state_identifier(df=macro_identify, var="deficit",countryname="Portugal", thresh=0.8, tol=0.01)
es1 <- state_identifier(df=macro_identify, var="deficit",countryname="Spain", thresh=0.8, tol=0.01)
se1 <- state_identifier(df=macro_identify, var="deficit",countryname="Sweden", thresh=0.8, tol=0.01)
uk1 <- state_identifier(df=macro_identify, var="deficit",countryname="United Kingdom", thresh=0.8, tol=0.01)
us1 <- state_identifier(df=macro_identify, var="deficit",countryname="USA", thresh=0.8, tol=0.01)

recdata2 <- as.data.frame(rbind(au1,at1,be1,ca1,dk1,fi1,fr1,de1,ie1,it1,jp1,pt1,es1,se1,uk1,us1))

recdata2 <- recdata2 %>% dplyr::rename(deficit_bin=state, deficit_weight=weight)

macro_data <- left_join(macro_data, recdata2, by=c("country","year"))

p1 <- ggplot(data=recdata2, aes(x=year, y=deficit_weight)) + geom_bar(stat = "identity", aes(y=deficit_bin), fill="grey", color=NA, alpha=1, width = 1.2) + geom_line(aes(x=year, y=deficit_weight), color="red", linetype=1, linewidth=1) + facet_wrap(facets="country") + xlab("Year") + ylab("Weight on Deficit Regime") + theme_pubclean() + theme(text = element_text(family = "Times New Roman"))

ggsave(p1, filename="deficit_regime.JPEG", height = 5, width = 9, dpi = 500)



###### differences across countries #########

macro_mean <- macro_lp %>% group_by(country,deficit_bin) %>%
  dplyr::mutate(deficit_mean = mean(deficit, na.rm=T)) %>% ungroup() %>%
  distinct(country,deficit_bin, .keep_all=T) %>% dplyr::select(country,deficit_bin,deficit) %>%
  dplyr::filter(!is.na(deficit_bin))

p1 <- ggplot(data=macro_mean, aes(x=country, y=deficit, fill=factor(deficit_bin))) + geom_bar(stat = "identity", position = "dodge") + labs(fill="Deficit Regime") + xlab("Country") + ylab("Average Budget Balance in % of GDP") + theme_pubclean() + theme(text = element_text(family="Times New Roman"), axis.line = element_line(color="grey", linewidth = 0.6), axis.ticks = element_line(color="grey"), legend.position = "top", axis.text.x = element_text(angle=50, hjust = 1.05, vjust = 0.9))

macro_mean <- macro_lp %>% group_by(country,recession) %>%
  dplyr::mutate(realgdpgr_mean = mean(realgdpgr, na.rm=T)) %>% ungroup() %>%
  distinct(country,recession, .keep_all=T) %>% dplyr::select(country,recession,realgdpgr) %>%
  dplyr::filter(!is.na(recession))

p2 <- ggplot(data=macro_mean, aes(x=country, y=realgdpgr, fill=factor(recession))) + geom_bar(stat = "identity", position = "dodge") + labs(fill="Recession") + xlab("Country") + ylab("Average GDP Growth Rate") + theme_pubclean() + theme(text = element_text(family="Times New Roman"), axis.line = element_line(color="grey", linewidth = 0.6), axis.ticks = element_line(color="grey"), legend.position = "top", axis.text.x = element_text(angle=50, hjust = 1.05, vjust = 0.9))

ggsave(p1, filename="deficit_mean_regime.JPEG", dpi=500, width=5, height=4)
ggsave(p2, filename="recession_mean_regime.JPEG", dpi=500, width=5, height=4)


####### robustness with 10% threshhold #########


## Real GDP Growth to identify slumps ##

au1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Australia", thresh=0.9, tol=0.01)
at1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Austria", thresh=0.9, tol=0.01)
be1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Belgium", thresh=0.9, tol=0.01)
ca1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Canada", thresh=0.9, tol=0.01)
dk1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Denmark", thresh=0.9, tol=0.01)
fi1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Finland", thresh=0.9, tol=0.01)
fr1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="France", thresh=0.9, tol=0.01)
de1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Germany", thresh=0.9, tol=0.01)
ie1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Ireland", thresh=0.9, tol=0.01)
it1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Italy", thresh=0.9, tol=0.01)
jp1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Japan", thresh=0.9, tol=0.01)
pt1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Portugal", thresh=0.9, tol=0.02) # sometimes have to increase tolerance
es1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Spain", thresh=0.9, tol=0.01)
se1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="Sweden", thresh=0.9, tol=0.01)
uk1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="United Kingdom", thresh=0.9, tol=0.01)
us1 <- state_identifier(df=macro_identify, var="realgdpgr",countryname="USA", thresh=0.9, tol=0.01)

recdata1 <- as.data.frame(rbind(au1,at1,be1,ca1,dk1,fi1,fr1,de1,ie1,it1,jp1,pt1,es1,se1,uk1,us1))

ggplot(data=recdata1, aes(x=year, y=weight)) + geom_line() + facet_wrap(facets="country")

recdata1 <- recdata1 %>% dplyr::rename(recession_10t=state) %>% dplyr::select(!weight)

macro_lp <- left_join(macro_lp, recdata1, by=c("country","year"))


## Deficits ##

au1 <- state_identifier(df=macro_identify, var="deficit",countryname="Australia", thresh=0.9, tol=0.01)
at1 <- state_identifier(df=macro_identify, var="deficit",countryname="Austria", thresh=0.9, tol=0.01)
be1 <- state_identifier(df=macro_identify, var="deficit",countryname="Belgium", thresh=0.9, tol=0.01)
ca1 <- state_identifier(df=macro_identify, var="deficit",countryname="Canada", thresh=0.9, tol=0.01)
dk1 <- state_identifier(df=macro_identify, var="deficit",countryname="Denmark", thresh=0.9, tol=0.01)
fi1 <- state_identifier(df=macro_identify, var="deficit",countryname="Finland", thresh=0.9, tol=0.01)
fr1 <- state_identifier(df=macro_identify, var="deficit",countryname="France", thresh=0.9, tol=0.01)
de1 <- state_identifier(df=macro_identify, var="deficit",countryname="Germany", thresh=0.9, tol=0.01)
ie1 <- state_identifier(df=macro_identify, var="deficit",countryname="Ireland", thresh=0.9, tol=0.01)
it1 <- state_identifier(df=macro_identify, var="deficit",countryname="Italy", thresh=0.9, tol=0.01)
jp1 <- state_identifier(df=macro_identify, var="deficit",countryname="Japan", thresh=0.9, tol=0.01)
pt1 <- state_identifier(df=macro_identify, var="deficit",countryname="Portugal", thresh=0.9, tol=0.01)
es1 <- state_identifier(df=macro_identify, var="deficit",countryname="Spain", thresh=0.9, tol=0.01)
se1 <- state_identifier(df=macro_identify, var="deficit",countryname="Sweden", thresh=0.9, tol=0.01)
uk1 <- state_identifier(df=macro_identify, var="deficit",countryname="United Kingdom", thresh=0.9, tol=0.01)
us1 <- state_identifier(df=macro_identify, var="deficit",countryname="USA", thresh=0.9, tol=0.01)

recdata2 <- as.data.frame(rbind(au1,at1,be1,ca1,dk1,fi1,fr1,de1,ie1,it1,jp1,pt1,es1,se1,uk1,us1))

recdata2 <- recdata2 %>% dplyr::rename(deficit_bin10t=state) %>% dplyr::select(!weight)

macro_lp <- left_join(macro_lp, recdata2, by=c("country","year"))


###### robustness LPs ########


lpb1 <- lp_builder(data = macro_lp, panelid = c("country","year"), endog = "socexp_public",
                   shock = "SPENDING", hor = 8, cum_mult = TRUE, vcov_mat = "vcovBK",
                   model = "within", effect = "individual", c_exog = c("realgdpgr","unemp"),
                   lag_exog = c("realgdpgr","unemp"), l_lags = 2, statedep = TRUE, statevar = "deficit_bin10t")

p1 <- plot_lp(lpb1, hor = 8, col1 = NULL, col2 = NULL, pointest = T, statedep = TRUE,
              statelabs = "High Deficits", ribbon = T, lty1 = NULL, lty2 = NULL,
              estimat = "within", lims = c(-2.7,0.75)) +
  theme(text = element_text(family="Times New Roman"), axis.line = element_line(color="grey", linewidth = 0.6), axis.ticks = element_line(color="grey"), legend.position = "bottom")

ggsave(plot = p1, filename = "IRF_DEFCRISIS10.JPEG", dpi = 800, height = 3, width = 3)


lpb2 <- lp_builder(data = macro_lp, panelid = c("country","year"), endog = "socexp_public",
                   shock = "SPENDING", hor = 8, cum_mult = TRUE, vcov_mat = "vcovBK",
                   model = "within", effect = "individual", c_exog = c("realgdpgr","unemp"),
                   lag_exog = c("realgdpgr","unemp"), l_lags = 2, statedep = TRUE, statevar = "recession_10t")

p2 <- plot_lp(lpb2, hor = 8, col1 = NULL, col2 = NULL, pointest = T, statedep = TRUE,
              statelabs = "Recession", ribbon = T, lty1 = NULL, lty2 = NULL, estimat = "within",
              lims = c(-2.6,0.75)) + theme(text = element_text(family="Times New Roman"),
                                           axis.line = element_line(color="grey", linewidth = 0.6),
                                           axis.ticks = element_line(color="grey"), legend.position = "bottom")

ggsave(plot = p2, filename = "IRF_RECESSION10.JPEG", dpi = 800, height = 3, width = 3)


##### cross-country identification ######

quantile(macro_lp$deficit, na.rm=T)

macro_lp$deficit_cc <- ifelse(macro_lp$deficit < -5, 1, 0)

cor(macro_lp$deficit_bin, macro_lp$deficit_cc, use = "na.or.complete")

quantile(macro_lp$realgdpgr, na.rm=T)

macro_lp$recession_cc <- ifelse(macro_lp$realgdpgr < 1, 1, 0)

cor(macro_lp$recession, macro_lp$recession_cc, use = "na.or.complete")

###### 

lpb1 <- lp_builder(data = macro_lp, panelid = c("country","year"), endog = "socexp_public",
                   shock = "SPENDING", hor = 8, cum_mult = TRUE, vcov_mat = "vcovBK",
                   model = "within", effect = "individual", c_exog = c("realgdpgr","unemp"),
                   lag_exog = c("realgdpgr","unemp"), l_lags = 2, statedep = TRUE, statevar = "deficit_cc")

p1 <- plot_lp(lpb1, hor = 8, col1 = NULL, col2 = NULL, pointest = T, statedep = TRUE,
              statelabs = "High Deficits", ribbon = T, lty1 = NULL, lty2 = NULL,
              estimat = "within", lims = c(-2.7,0.75)) +
  theme(text = element_text(family="Times New Roman"), axis.line = element_line(color="grey", linewidth = 0.6), axis.ticks = element_line(color="grey"), legend.position = "bottom")

ggsave(plot = p1, filename = "IRF_DEFCRISIS_CC.JPEG", dpi = 800, height = 3, width = 3)


lpb2 <- lp_builder(data = macro_lp, panelid = c("country","year"), endog = "socexp_public",
                   shock = "SPENDING", hor = 8, cum_mult = TRUE, vcov_mat = "vcovBK", 
                   model = "within", effect = "individual", c_exog = c("realgdpgr","unemp"),
                   lag_exog = c("realgdpgr","unemp"), l_lags = 2, statedep = TRUE, statevar = "recession_cc")

p2 <- plot_lp(lpb2, hor = 8, col1 = NULL, col2 = NULL, pointest = T, statedep = TRUE,
              statelabs = "Recession", ribbon = T, lty1 = NULL, lty2 = NULL, estimat = "within",
              lims = c(-2.6,0.75)) +
  theme(text = element_text(family="Times New Roman"), axis.line = element_line(color="grey", linewidth = 0.6), axis.ticks = element_line(color="grey"), legend.position = "bottom")

ggsave(plot = p2, filename = "IRF_RECESSION_CC.JPEG", dpi = 800, height = 3, width = 3)


#######

write_dta(macro_data, "as_macro_data.dta")
write_sav(macro_data, "as_macro_data.sav")
saveRDS(macro_data, "as_macro_data.RData")





